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SUMMARY 

The next generation of commercial aircraft will include 
turbofan engines with performance levels significantly better 
than those in the current fleet. Control of particulate and 
gaseous emissions will also be an integral part of the engine 
design criteria. These performance and emission requirements 
present a technical challenge for the combustor: control of the 
fuel and air mixing and control of the local stoichiometry will 
have to be maintained much more rigorously than with com- 
bustors in current production. A better understanding of the 
flow physics of liquid fuel spray combustion is necessary. 
This paper describes recent experiments on spray combustion 
where detailed measurements of the spray characteristics were 
made, including local drop-size distributions and velocities. 
Also, an advanced combustor CFD code has been under deve- 
lopment and predictions from this code are compared with 
experimental results. Studies such as these will provide 
infonnalion to the advanced combustor designer on fuel spray 
quality and mixing effectiveness. Validation of new fast, 
robust, and efficient CFD codes will also enable the combus- 
tor designer to use them as valuable additional design tools for 
optimization of combustor concepts for the next generation of 
aircraft engines. 

A. INTRODUCTION 

Aircraft engines being envisioned now for the next-century 
aircraft will have requirements that present formidable techni- 
cal challenges to the combustor designer. In the subsonic 
commercial transport arena, demand for low operating cost 
translates into reduced fuel consumption and improved dura- 
bility and reliability. Higher operating pressures of the 
combustor are forecast, with higher resulting fuel system 
turndown ratios. Coupled with these performance requirements 
is the demand world wide for control of pollutant emissions 
from aircraft engines, especially oxides of nitrogen. In the 
arena of supersonic commercial transports, forecasts are 


predicting fleets of hundreds of transports operating at 
Mach 2.0 to 2.5, with airfares only slightly higher than today's 
long-range subsonic fares. Economical fuel consumption is a 
requirement. But the largest technical challenge here is low 
levels of oxides of nitrogen emissions during engine cruise 
such that there would be no adverse impact to the earth's 
environment, specifically the atmosphere's ozone layer. 

These technical and environmental challenges represent design 
requirements for the combustor that are outside of engine 
companies' experience. Empirically based design methods are 
insufficient by themselves. To augment this design system, the 
companies are increasingly turning to computational fluid 
dynamics (CFD) computer codes. Severe limitations with cur- 
rently available codes are excessive time requirements to run 
a CFD code to analyze a complex combustor design: a penalty 
both in turn-around time for design answers as well as cost for 
the calculations. We at NASA are attempting to improve this 
current situation in the industry by developing a fast, robust, 
efficient computer code for internal chemical reacting flows. 
The objective is to produce a CFD code that can be used as 
a more powerful design tool in the industry to analyze com- 
plex combustor designs in significantly shorter turn-around 
time than current computer codes can achieve. 

An integral part of the development of a new CFD code is the 
validation of this code with experiments that represent the 
complex features of the flows which need to be analyzed. 
Therefore, at NASA we are also conducting experiments on 
liquid fuel spray combusting flows with increasingly complex 
features. This will provide some of the required data for 
validation of our CFD code and other's codes. Detailed mea- 
surements of both liquid sprays and gas characteristics are 
being obtained under both nonburning and burning conditions. 

This paper will highlight recent results that NASA has 
obtained with its spray combustion experiment and describe 
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the formulation and physical modeling of its new spray com- 
bustion CFD code (ALLSPD). The application of this code to 
combustion problems will be illustrated by several examples. 

B. SPRAY COMBUSTION EXPERIMENTS 

Combusting sprays are very important for gas turbine engine 
applications. The investigation of combusting spiays should 
lead to a better understanding of die physics involved in this 
complicated process. Important processes involved in combust- 
ing sprays are the interactions between the droplets and the 
gas phase, the vaporization of the droplets, and chemical 
reactions with heat release. These physical processes are 
coupled and can only be completely described using numerical 
modeling. As part of an effort to improve the numerical 
modeling of gas turbine combustors, an experimental study 
has been undertaken to obtain a data set for a relatively simple 
liquid-fueled combustor that can be used for comparison with 
numerical models (Ref. 1). 

Because of their practical applications, swirling flows with 
combustion have been studied by a large number of investi- 
gators. Earlier reviews of swirling flows both with and without 
combustion present some general trends (Refs. 2 to 4). These 
papers predate die development of nonintrusive. laser-based 
diagnostics; consequently all of the results described were 
obtained using intrusive instrumentation and detailed structure 
measurements for these types of flows were not possible. With 
the advent of newer instrumentation techniques, namely laser 
Doppler anemometry. additional details of die structure of 
these types of flows began to emerge. Laser Doppler anemo- 
metry velocity measurements in spray flames (Refs. 5 and 6) 
reveal some of the flowfield structure of swirling flames. The 
development of the phase/Doppler particle analyzer (Ref. 7). 
enabled the simultaneous measurement of droplet size and 
velocity. This instrument has been used by a number of inves- 
tigators for measurements in spray flames in a variety of 
configurations (Refs. 8 to 12). This instrument has the 
capability to measure velocities of both the gas and droplet 
phases in a combusting spray. 

B.l EXPERIMENT 

The combustor utilized in the present experiment is illustrated 
in Fig. I. It consists of a center mounted air-assist fuel nozzle, 
Parker Hannifin research simplex air-assist atomizer, sur- 
rounded by a coflowing air stream. The nozzle orifice diame- 
ter was 4.8 mm. Both the air assist and the coflow air streams 
had swirl imparted to them using 45 degree swirlers. The 
swirlers were constructed by machining 45 degree slots into 
rings. Both streams were swirled in the same direction for the 
present study. The combustion air was not preheated and 
entered the combustor at 297 K. The top of the air-assist 
nozzle was water cooled to prevent overheating of an O-ring 
in the nozzle assembly. The temperatures of the fuel, atomiz- 
ing air and coflow air streams were measured using Giromal 
Alumel thermocouples. Flow rates of the air streams were 
measured using calibrated orifices and the fuel flow rate was 
measured using a mass flowmeter. All results reported in the 
present study are reported for a coflow air flow rate of 
13.88 g/s, an air-assist flow rate of 0.96 g/s. and a fuel flow 
rate of 0.38 g/s. The fuel used was heptane. The coflow 


stream entered the combustor in three radial locations, passed 
through a honeycomb flow straightener, and the swirlers 
before exiting the combustor. The swirler was located 140 mm 
upstream of the combustor exit. The flow from the combustor 
discharged into ambient, stagnant surroundings. 

The combustor was mounted vertically within a large 
(1.8 m sq by 2.4 m high) enclosure. The entire enclosure was 
mounted on two sets of linear bearings and was traversed 
using stepper motors to provide motion in two directions. The 
combustor assembly itself could be traversed in the vertical 
direction using a third stepper motor to allow measurements 
at all locations in the flowfield. This arrangement allowed 
rigid mounting of all optical components. 

The phase/Doppler particle analyzer was used for all mea- 
surements reported in this study. A schematic of the two- 
component instrument is shown in Fig. 2. The beam from a 
6 W Argon-Ion laser is split into 488.0 and 514.5 mn wave- 
lengths using a dichroic mirror. Each beam is then focused 
onto a rotating diffraction grating which splits each beam into 
several pairs. The two first-order beams for each wavelength 
are then recombined onto the optical axis using a dichroic 
mirror, collimated and focused at a point to form the two- 
component probe volume. In the present study, the transmit- 
ting optics utilized a 500 mm focal length lens. The receiving 
optics were located 30 degrees off axis in the forward-scatter 
direction. Light was collected using a 500 mm focal-length 
lens and then focused onto a 100 pm by 1 mm long slit. The 
collected light is then split and picked up by four photo- 
detectors. Three are arranged to look at the signals from the 
514.5 beams and one receives light from the 488 nm beams. 
Each of the three photodetectors for the green beams are 
imaged at a different area of the collection lens and the phase 
difference between the signals is used for the size determina- 
tion. Details of the instrument can be found in Ref. 7. 

In the present study, velocities of both the liquid and gaseous 
phases were measured. This was accomplished by seeding the 
gas phase with nominal 1 pm size aluminum-oxide particles. 
The coflow, air-assist flow, and the ambient surroundings 
were all seeded to minimize biasing. Phase discrimination is 
inherent in the instrumentation with the ability to size each 
measured particle. At each spatial location, two measurements 
were taken in order to accurately measure the velocity of each 
phase. A threshold voltage for die photodetectors at the speci- 
fied laser power was determined experimentally, below which 
signals from die aluminum-oxide particles were not detected. 
For the droplet measurements, the photodetector voltage was 
kept below this threshold value in order to eliminate interfer- 
ence from the aluminum-oxide particles. Total laser power for 
all wavelengths was fixed at 1 .5 W for all the measurements. 
Particles with diameters less dian 2.4 pm were used to repre- 
sent the gas phase velocity. Two complete traverses were 
taken in order to measure all three components of velocity and 
provide a check on flow symmetry. Each traverse measured 
axial velocity and either radial or angular velocity. Generally. 
64 000 measurement attempts were made at each measurement 
location. The percentage of measurements actually validated 
depended on the number density and velocities of drops at 
each location and ranged from about 65 to 90 percent. 
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B.2 RESULTS AMD DISCUSSION 

In the present study, results are presented for a single axial 
location at 5 inni downstream of the nozzle. Gas phase results 
for mean velocities are presented in Figs. 3(a) to (c) for 
isothermal, single-phase flow without droplets and two-phase 
flow with combustion. Mean gas phase axial velocity, 
presented in Fig. 3(a). presents results from a complete 
traverse across the combustor and illustrates the symmetry of 
the flowfield. The combustor exit dimensions are illustrated on 
the x-axis of the figure for reference. A small recirculation 
zone is evident near the center of the nozzle. At this axial 
location relatively close to the nozzle, velocity gradients are 
extremely large in the flow from the air-assist stream contain- 
ing the droplets. Effects of combustion on the flowfield are 
significant. Both the maximun and minimum mean axial 
velocities are increased for the combusting case compared to 
the isothermal case due to the gas expansion associated with 
the heat release. For the case with combustion, velocities 
increase from nearly 0 to 3.*) m/s and then decrease to -30 m/s 
within a radius of about 12 mm. Mean axial velocities in the 
coflow stream are not affected by the combustion at this axial 
location. 

Figure 3b presents mean radial velocity for the gas phase. As 
shown in Fig. 3(b). effects of combustion are very dramatic 
for radial velocity. Maximum radial velocities increased ffom 
about 1 0 m/s for the isothermal case to about 40 in/s for the 
case with combustion due to the radial expansion of the gas. 
Again, the gas from the coflow stream is not affected by the 
combustion at (his axial location. 

Mean gas phase angular velocities are presented in Fig. 3(c). 
For this case, reaction and the presence of droplets decreases 
the maximum angular velocities in the flowfield. Some of the 
decrease in angular velocity for the gas phase can be attri- 
buted to the momentum transferred to the droplets since they 
do not initially have a swirl component. 

Fluctuating gas phase velocities are presented in Figs. 4(a) to 
(c) for both the single-phase, isothermal and the two-phase, 
combusting cases. All fluctuating velocities presented are root- 
mean-squared (rms) values. Figure 4(a) presents radial profiles 
of gas phase fluctuating axial velocity. Maximum values of 
fluctuating axial velocity are similar for both the combusting 
and isothermal cases. The case with combustion does show 
larger values of axial rms velocity at radial locations between 
approximately 5 and 15 mm from the center of the nozzle. 
Axial velocities are also higher at these locations for the 
combusting case, see Fig. 3(a). Fluctuating radial velocities, 
illustrated in Fig. 4(b). show dramatic differences between the 
isothermal and combusting cases. The maximum velocity loca- 
tions have shifted radially outward corresponding to the shift 
in mean radial velocity, see Fig. 3(b). The maximum fluctuat- 
ing radial velocity has also increased from about 10 nv's to 
15 nv's. 

Fluctuating angular gas phase velocities are presented in 
Fig. 4(c). Similar to the results shown for mean angular gas 
phase velocities, fluctuating angular velocities generally 
decreased with combustion and the presence of the liquid 
phase compared to Die single-phase, isothermal case. A small 
region from a radius of about 7 to 15 mm shows increased 
values of fluctuating angular velocity for the case with combustion. 


Mean velocities for the drops are presented in Figs. 5(a) to (c) 
for the case with combustion. In the experimental study, 
velocities were measured for drop sizes ranging ffom 4 to 
142 pm. Results are presented for drop sizes of 15. 32 and 
52 pm. Measured gas phase velocities are also presented in 
the figure. Note that results are only illustrated from -15 to 
+ 1 5 mm for the radial direction because no drops were present 
at larger radial locations. Figure 5(a) presents mean drop axial 
velocity at 5 min downstream. Similar to the results previously 
shown for Die gas phase, the flowfield is very symmetric. 
Axial velocity is correlated with drop size in all regions. In 
the main region of the spray, at a radius of about 7 mm, the 
maximum velocity of the gas phase was about 38 nv's. and 
about 28 in/s for the 32 pm drops. Even the maximum veloc- 
ity of the 15 pm drops lagged the gas phase by about 5 m/s. 
In the center of the flowfield is a small recirculation zone, see 
Fig. 3(a). There, only the 15 pm drops showed negative axial 
velocities while larger drops had positive velocities. 

Mean drop radial velocities are presented in Fig. 5(b). Again, 
there is a correlation between drop size and velocity in the 
flowfield. Maximum mean radial velocities are slightly higher 
than maximum axial velocities for the drops due to the heat 
release and radial expansion of the gas. Mean angular veloci- 
ties of the drops are presented in Fig. 5(c). Angular velocity 
is not as symmetric and is also much smaller than the other 
two components of velocity. The mean drop angular velocity 
is a strong function of the drop size with the smaller drops 
showing the least velocity difference with the gas phase. 

Fluctuating droplet axial, radial, and angular velocity com- 
ponents for the three drop sizes and gas phase are presented 
in Figs. 6(a) to (c). respectively. The fluctuating drop veloci- 
ties presented are root-mean-squared (rms) values. Generally, 
the smaller drops are affected more by the gas phase turbu- 
lence and have larger fluctuating velocities than the larger 
drops. Velocity fluctuations are clearly not isotropic since 
fluctuating axial and radial velocities are considerably larger 
than fluctuating angular velocities. 

In addition to drop mean and fluctuating velocities, the liquid 
volume flux is important in two-phase flows. Drop number- 
flux measurements are presented in Fig. 7, where results are 
illustrated for four drop size groups. As shown in Fig. 7, 
number flux is nearly symmetric. The results show that 
smaller droplets have much larger number fluxes. The distri- 
bution of the larger droplets is still very important since much 
of the liquid mass is contained in the larger droplets. 
Relatively few drops are found in the center region of the 
flowfield due to Die 45 degree swirler that is used in the air- 
assist stream. 

C. COMBUSTION CFD CODE 

The objective of the present work is to develop a numerical 
solution procedure which can efficiently handle the coupling 
between a spray model and a well-developed strongly implicit 
flow solution algorithm (Refs. 13 and 14). In Die past, spray 
models have been coupled with different flow algorithms and 
shown some promising results (Refs. 15 to 17). However, 
most of the spray models were coupled with a flow solver 
employing a segregated approach, such as a TEACH-type 
code (Refs. 15 and 16). which has been used very extensively 
in the industries for the past two decades. Although simple 
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and easy to implement, the TEACH-type code usually suffers 
poor convergence due to the explicit (or semi-implicit) 
treatments of the chemical source terms and the sequential 
solution approach. Recent development of CFD techniques and 
the advent of computer technology have allowed us to explore 
more ambitious schemes to solve reacting flow problems. 
Strongly coupled and implicit numerical schemes, although 
requiring much more computer storage and complexity of the 
algorithm, have been very popular for nonreacting compressi- 
ble (Refs. 18 and 19) and incompressible (Refs. 20 and 21) 
flow computations. 

Shuen and Yoon (Ref. 22) developed a coupled scheme for 
high speed reacting flows. RPLUS. which has been used and 
studied quite extensively in recent years. However, like other 
compressible flow solution algorithms. RPLUS may not be 
suitable for low speed flow computations. There are two well- 
recognized reasons (Refs. 23 to 25) for the convergence 
difficulties related to compressible flow codes. First, the 
system's eigenvalues become stiff at low flow velocities. 
Second, the pressure term in the momentum equation becomes 
singular as the Mach number approaches zero, which yields a 
large roundoff error and smears the pressure variation field. 
This will not only result in slow convergence but often will 
produce inaccurate solutions. To circumvent Hie above diffi- 
culties. Shuen et al. (Ref. 13) developed a coupled numerical 
algorithm for chemical nonequilibrium viscous flows. 
ALLSPD. which utilizes Hie decomposition of the pressure 
variable into a constant reference pressure and a gauge 
pressure to reduce the roundoff errors and adds a precondition- 
ing time derivative term to rescale the system eigenvalues. 
The results of Hiese treatments show that the convergence 
properties are almost independent of the flow Mach number. 

Numerous spray models have been proposed and investigated 
for different spray combustion problems (Refs. 26 and 27) in 
the past decade. Recent spray models differ in specific details, 
but generally may be divided into two categories: locally 
homogeneous flow (LHF) models and separated flow (SF) 
models. LHF models represent the simplest treatment of a 
multiphase flow and have been widely used to analyze sprays 
(Ref. 28). The key assumption of the LHF model is that inter- 
phase transport rates are fast in comparison to the rate of 
development of the flow. This implies that all phases have 
identical properties at each point in the flow. Clearly. LHF 
models are only formally correct for flows containing infi- 
nitely small droplets. 

Numerous SF models have been proposed to consider inter- 
phase transport phenomena (Ref. 29). Among diem, the dis- 
crete droplet approach (Refs. 26 and 30) has been adopted, 
since it reduces numerical diffusion while providing a conve- 
nient framework for dealing with multiple droplet size and 
complex interphase transport phenomena. Many discrete drop- 
let models neglect the effects of turbulence on interphase 
transport (Refs. 26 and 30). This implies that droplets follow 
deterministic trajectories, yielding the deterministic separated 
flow (DSF) model. Neglecting the effects of turbulence on 
droplet transport is appropriate when characteristic droplet 
relaxation times are large in comparison to characteristic times 
of turbulent fluctuations. Few practical sprays, however, 
satisfy this condition. Dukowicz (Ref. 31) and Gosman and 
loannides (Ref. 32) have adopted stochastic methods to study 
droplet dispersion by turbulence. Faeth and coworkers 


(Ref. 33) extended the analysis of Gosman et al. to include the 
effects of turbulence on interphase heat and mass transport. 
Their stochastic separated flow (SSF) model (Ref. 33) has 
been evaluated in a wide variety of parabolic flows with very 
encouraging results. 

In a recent paper (Ref. 34), Hie ALLSPD algorithm was 
extended to include a SSF spray model (Ref. 33), a recently 
developed low Reynolds number K-e turbulence model 
(Ref. 35) and a multiblock treatment to calculate the gas 
turbine combustion flows, where the liquid spray is an im- 
portant ingredient of combustion. Although the turbulence 
model has been included in this study for turbulent combus- 
tion flow calculations, the mean flow quantities are still used 
in Hie chemistry calculations. The turbulence closure problem 
for chemistry is more complicated and computationally inten- 
sive. The consideration of a suitable turbulent combustion 
closure model is left for our next-phase study in the near 
future. 


In this paper, in addition to the gas turbine spray combustion 
flow, a single-phase premixed turbulent combustion flow and 
a nonreacting turbulent flow are also included to demonstrate 
the current status of this research. In the following sections, a 
brief mathematical formulation of the governing equations for 
both gas and liquid-phases is described. The numerical method 
and discretization procedure are given next and finally some 
sample numerical results are presented. 


C.l GOVERNING EQUATIONS 

C.la Gas-Phase Equations 
C.laO) Navier-Stokes Formulation 

The two-dimensional, unsteady, compressible, density- 
weighted time-averaged Navier-Stokes equations and species 
transport equations for a chemically reacting gas of N species 
written in generalized nonorthogonal coordinates can be 
expressed as 
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In the above expressions, r, 5, and q are the time and spatial 
coordinates in the generalized coordinates and E, and q, are 
the grid speed terms. The E v . q x , and q v are Hie metric 
terms and the J is Hie transformation Jacobian. The power. 6. 
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is an index for two types of governing equations with 5 = 0 
for two-dimensional and 5 = I for axisymmelric cases (with 
x being the axial and v the radial coordinates, respectively). 
The vectors Q, E. F, F^, F v , H c and H ( in the above defini- 
tions are 
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where p. p. it, v. K, c. and Y/ represent the density, pressure. 
Cartesian velocity components, turbulent kinetic energy, 
dissipation rate of turbulent kinetic energy and species mass 
fraction, respectively; E = e + H(u 2 + v 2 ) is the total 
internal energy with e being the thermodynamic internal 
energy; and S, is the rate of change of species i due to 
chemical reactions. The normal and shear stresses, energy, 
species, and turbulent diffusion fluxes are given by 


, du 2 

S ~ 3 


T .tv “ 


du dv 

izy 

- dv 2 
2M < W ' 3 M ‘ 


du + dv' 


du + dv' 
dx 


4 v 2 

Tf * ^ 3 7 ‘ 3 


1 du dv ' 

W ~ * 

* / 

ar, 


■ *' -ST * p S " lD "’ TT 

q " ’ *• Sp * p S 17 






M, 

M, + 

o. 


0K 

"ST 

0K 

5v 

dc 

aT 


6 


T 


ye 




0e 

~dy 


dY i 

</.* - pD, m . 


dY l 

Q\1 m P • 


where T, pi^, p r pi,. and k e are die temperature, effective vis- 
cosity. turbulent viscosity, molecular viscosity and the effec- 
tive thermal conductivity, respectively; 
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is the effective binary diffusivity of species i in tlie gas 
mixture. X, the molar fraction of species i. and D t j the binary 
mass diffusivity between species / and j. The quantities related 
to the source term in the turbulent equations are given as 
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where y„ in tlie expression of R f is the normal distance away 
from the wall, k t is the molecular thermal conductivity, C pm 


is the specific heat of the gas mixture and Pr, is the turbulent 
Prandtl number. 


The vector H, in Eq. (I) represents the source term that 
accounts for the interactions between the gas and liquid 
phases. In vector H,. n is the number of droplets in the /;th 
characteristic group of droplets; m is the evaporation rale 
of each particle group; p p is the liquid density; r p is the 
droplet radius; n p and v p are the particle velocities at the /;lh 
group; lip, is the enthalpy of fuel vapor at the droplet surface 
and I, AT is the convective heat transfer between two phases. 
Detailed discussions about the liquid phase equations will be 
described later. The temperature and pressure are calculated 
iteratively from tlie following equations 

c = £ Yfr - L. I,,- I,; * J r r CftdT. 


* p - p R J E 

i-i w t 

where R„ and T yc j are the universal gas constant and reference 
temperature for thermodynamic properties, and IV,. r,. /i,.//^ 
are the molecular weight, constant pressure specific heat, 
thermodynamic enthalpy, and heat of formation of species /, 
respectively. 

In reacting flow calculations, the evaluation of thermophysical 
properties is of vital importance. In this papier, the values of 
C pl , k tl , and p„ for each species are determined by fourth- 
order polynomials of temperature, as described in Shuen 
(Ref. 36). The specific heat of the gas mixture is obtained by 
mass concentration weighting of individual spiecies. The 
thermal conductivity and viscosity of the mixture, however, 
are calculated using Wilke's mixing rule (Ref. 37). The binary 
mass diffusivity £>„ between species i and j is obtained using 
the Chapman-Enskog theory (Ref. 37). 


C.la(2) All-Mach-Number Formulation 

As noted earlier in the introduction section, the two main 
difficulties tliat render tlie compressible flow algorithms 
ineffective at low Mach numbers are the roundoff error caused 
by die singular pressure gradient term in the momentum equa- 
tions (the pressure term is of order 1/M 2 while the convective 
term is of order unity in tlie nondimensional momentum equa- 
tions) and the stiffness caused by the wide disparities in 
eigenvalues. To circumvent these two problems regarding the 
low Mach number calculations, following tlie approach by 
Merkie and Choi (Ref. 25), Shuen et al. (Ref. 13) added a 
time preconditioning term to rescale the system eigenvalues 
and decompose the pressure variable into a constant reference 
pressure part and a gauge pressure part. This all-Mach-number 
formulation has been extended to include the turbulent and 
spray equations. The resulting Navier-Stokes equations in a 
conservative form are 



5Q + f (E - E v ) + 0(F - K) 
~5T — 35 — — 3n — 


U c * H,. 

(3) 


where the primitive variable vectorQand the preconditioned 
matrix f are given as 
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where T* is the pseudotime. (3 is a parameter for rescaling the 
eigenvalues of the new system of equations: 
and/t = e + p / pis the specific enthalpy of die gas mixture. 
The definition of vectors in Eq. (3) is identical to those in 
Eq. (1) except that the pressure terms in the momentum 
equations are replaced by gauge pressure. p g . The derivation 
of the al!-Mach-number formulation can be found in detail in 
Shuen et al (Refs. 13 and 14). 


C.lb. Liquid-Phnse Equations 
C.lb(I) Droplet Motion Equations 


The liquid phase is treated by solving Lagrangian equations of 
motion and transport for the life histories of a statistically 
significant sample of individual droplets. This involves 
dividing the droplets into n groups (defined by position, 
velocity, temperature and diameter) at the fuel nozzle exit and 
then computing their subsequent trajectories in the flow. The 
spray model used in this study is based on a dilute spray 
assumption which is valid in the regions of spray where the 
droplet loading is low (Refs. 17 and 33). Tiie liquid fuel is 
assumed to enter the combustor as a fully atomized spray 
comprised of spherical droplets. The present model does not 
account for tiie effects due to droplet breakup and coalescence 
processes which might be significant in a dense spray situa- 
tion. The Lagrangian equations governing the droplet motion 
are 


d.x 

dr r 

(4) 

<h' r 

■ v_ , 

(5) 

(It r 


. 3 W'Vl, 
lit 16 0 s l * PI 

v n p 

(6) 

dv P „ 3 C D p f Re p / _ , 

dt 16 i pl ' 

(7) 

pp'p 


where the particle Reynolds number, Re p . and 
coefficient. Cp. are defined as 

the drag 


Re P a 2 ~ (("* ■ “pf + h ■ v rf 


t r. 


24 


R ‘n 

0.44 


1 + 


Re 


2/3 


for Re p < 1000. 


for Re,, > 1000. 


The subscript g represents the gas-phase quantities and p 
represents the liquid-phase (or 'particle") quantities. Equations 
(4) and (5) are used to calculate the new droplet positions and 
Eqs. (6) and (7) are used to update tiie new droplet velocities 
at the new droplet locations. A second-order Runge-Kutta 
scheme was used to integrate Eqs. (4) to (7). 


C.lb(2) Droplet Mass and lieat Transfer Equations 

As described by Faeth (Ref. 26), the following correlations 
were used to approximate the mass and heat transfer coeffi- 
cients for a single isolated spherical droplet: 


th„ <L , . 

r r - 2AUn I + B), 
P D f 


( 8 ) 


hd p _ 2N r ln(l + B) u ~' 

— 1(1 * sf' - ij ' 


(9) 


where /; is the heat transfer coefficient, d is the droplet 
diameter. is the fuel mass diffusivity. k is the thermal 
conductivity of fuel vapor and rti ” is the fuel mass evapora- 
tion rate per unit area. The N x ana N are defined as 


N . 


1 + 


027 6 Re ^ Pi 1/3 
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1.232 
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OnORe'^Sc 1/3 


1 + 


1.232 
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Re p Sc 

where 5c and Le are the Schmidt and Lewis numbers, respec- 
tively. The Spalding number, B, in Eqs. (8) and (9) is defined 




Y fw 


7* 


1 


( 10 ) 


where Y, vp is the fuel vapor mass fraction at the surface of the 
droplet and is the mean fuel mass fraction of the ambient 
gas. In the present study. Y ftp is obtained from the following 
equation 


l fgp 


Vr 


vF\' -*/»>•' 


(M) 
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where )V a is the molecular weight of gas excluding fuel vapor. 
W'f is the molecular weight of fuel and Xf gp is the mole 
fraction of fuel. The X^, is obtained from the assumption of 
Raoult's law. Based on this assumption, the mole fraction at 
the droplet surface is equal to the ratio of die partial pressure 
of fuel vapor to the total pressure. For the present spray 
calculation, the partial pressure of fuel vapor was computed 
based on the following empirical correlation (Ref. 37): 

t \ 
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K, 
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Py jX 
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+ P v yr + P v4 x' 


( 12 ) 


where x « 1 - T/7 ( „ P vl = -7.28936. P v2 - 1.53679. P v3 
- -3.08367. P v4 » -1.02456. T c and p c are the critical 
temperature and critical pressure of the fuel vapor, 
respectively. 


where the value of a is between 0 and 1 with a ■ 0 referring 
to the vortex center and ct-1 referring to the droplet surface. 
The initial and boundary conditions for Eq. (13) are 


'in]' 


‘in]' 


a “ 0. 
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where 1 

^ J, 

droplet surface 


is obtained from die energy balance at die 
by the following equation 


C.lb(3) Dropfet Internal Temperature Equations 


As a single droplet enters a hot environment, the immediate 
small portion of the droplet near the surface will be heating up 
quickly while the center core of the droplet remains "cold". 
The heat will be conducted and convected to the entire interior 
as the droplet penetrates further into the hot ambient gas. 
Eventually, the temperature within the droplet will become 
nearly uniform before the end of its lifetime. To solve this 
transient phenomena within the droplet is not trivial. In the 
past, certain approximations (Ref. 38) are usually made to 
alleviate this computational burden while obtaining reasonably 
good results. Among them, the simplest one is die uniform 
temperature model. This model assumes that the thermal con- 
ductivity of the fuel is infinite. Of course, this is not valid at 
the beginning of the heating-up process of the droplet. 
Another model considers the heat diffusion inside the droplet. 
This yields a conduction model. The temperature distribution 
within the droplet is obtained by solving the one-dimensional 
heat conduction equation subject to the convective boundary 
conditions at the droplet surface. The conduction model 
completely neglects the convective phenomena within the 
droplet which might occur due to significant shear forces at 
the surface induced by high slip velocities. Tong and 
Sirignano (Ref. 39) developed a vortex model for the internal 
temperature of a single droplet which accounts for the 
convective effect of the Hill's vortex formation inside the 
droplet. In this study, Tong and Sirignano's vortex model is 
applied to obtain the internal temperature distribution of the 
droplet. The equation governing the internal temperature 
distribution based on this model is 


dT r _ */ 


ST 


17 


c r^r 


d-T. , , ar 

a — p + (l + C(r) a) p 
da - da 


(13) 


and 


c«, - _ 


c r i p l 


dt 





AT - at p lift I 


(14) 


where rh p and h are calculated from Eqs. (8) and (9). respec- 
tively. hf g is the latent heat of the fuel and AT = 7^ - T pr . 
where T , f is droplet surface temperature and 7^, is the mean gas 
temperature evaluated in the following way 


j = 1 t + 2 t 

* I s J ps 


An implicit scheme was used to solve Eq. (13) subject to the 
initial and boundary conditions. Second-order central differ- 
ences were used for the spatial differential terms and a first- 
order difference for the time term. These treatments rendered 
a scalar tridiagonal algebraic system and was solved by the 
Thomas algorithm (Ref. 40). 


C.2 NUMERICAL METHODS 
C.2a Discretized Equations 

Equation (3) is the final gas-phase governing equation to be 
solved numerically. To obtain time-accurate solutions for 
time-evolving problems, a dual time-stepping integration 
method can be applied to Eq. (3). The solution converged in 
pseudotime corresponds to a time-accurate solution in physical 
time (Ref. 13). However, for the present study, since only the 
steady state solution is of our interest, the physical time term 
in Eq. (3) can be dropped and the solution can be marched 
completely in pseudotime to obtain the final steady state 
solution. One advantage of marching to the steady state 
solution in pseudotime is that the convergence of the marching 
(iterative) process is determined by the eigenvalue characteris- 
tics on the pseudotiine space and not by the original stiff 
eigenvalues. The analysis of the eigensystem has been per- 
formed m our previous study (Refs. 13 and 14) and therefore 
will not be repeated here. It should be noted that the inclusion 
of the turbulent x-c equations in this study does not affect 
the system eigenvalues at all and, therefore, the properties of 
the all-Mach-number formulation analysis in previous papers 
(Refs. 13 and 14) are still valid. 
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After linearization and applying a first-order time differencing, 
Eq. (3) can be expressed in the following form: 


T - At’D + Ar’ 


aA 


a 


HZ dZ 


+ At* 


dB 

h*n 


nn 3n 


* aO - - At*(R f. 


(15) 


Since a multiblock treatment is applied to the present 
numerical algorithm, the interfaces between blocks (or zones) 
become a special kind of boundary. At the interface, the 
governing equations are still solved without any distinction 
from the rest of the interior points. The flow variables have 
been carefully arranged in such a way that the information 
from the neighboring block(s) is automatically brought into 
the calculation at the interface of the current block. The 
multiblock treatment of the present study not only provides 
flexibility for complex geometry calculations but also reduces 
the size of tire storage array for tire MSIP coefficients due to 
the smaller sizes of the subdomains. 


where 


R r 


. a (* - ^ 


ar 


dn 




(16) 


where p denotes the previous iteration level. D is the Jacobian 
for chemical and turbulent source terms, A and B are the 
inviscid term Jacobians and R^ and R nn are the viscous term 
Jacobians. The expressions for these Jacobians can be found 
in Shuen et al. (Refs. 13 and 14) except for the turbulent part. 
Central differences were used to discretize the spatial deriva- 
tive terms in Eqs. (15) and ( 16) for both explicit and implicit 
operators. The resulting coupled algebraic equations are solved 
using a modified strongly implicit procedure (MSIP) proposed 
by Schneider and Zedan (Ref. 4 1 ). which is completely vector- 
izable along i + 2/ diagonal direction. 

C.2b Boundary Conditions 

The boundary conditions for gas phase equations are described 
as follows. Here, only the subsonic flow boundary conditions 
are considered in the present paper. At the inlet, all quantities 
are specified except pressure which is obtained through extra- 
polation from the pressure at interior points. At the exit, the 
governing equations are solved at the exit station by applying 
backward differences for die streamwise derivative terms 
(central differences are still used for cross-stream derivative 
terms). However, the streamwise pressure derivative terms are 
centrally-differenced. This treatment requires the pressure 
information at the station one step downstream (outside the 
computational domain) where a constant pressure condition is 
enforced. At the symmetry line, the governing equations are 
solved via the use of the symmetry conditions for two- 
dimensional flows. For die axisymmetrie case, a singularity 
exists at this line and a simple one-sided difference is used to 
implement the symmetry conditions. At the solid wall, no slip 
conditions are used for the velocities and an adiabatic wall is 
assumed. Normal derivatives for species. Y f . and (p + 2/3 pK) 
are set to zero. For the turbulent quantities, the low Reynolds 
number turbulence model (Ref. 35) used in this study requires 
the specification of the k and c at the walls as follows: 

k = 0.250 M ' 

4 

“t 

e = 0.251 _L 

v 

where u T is the friction velocity at the wall. 


C.2c Spray Source Terms 

In Eq. (3), the liquid-phase interaction with the Navier-Stokes 
equations is modeled as a source term which accounts for 
mass, momentum and energy exchanges between two phases. 
In order to compute this source term, the liquid-phase govern- 
ing equations described in the previous section are integrated 
in time once initial conditions have been specified. As a liquid 
droplet begins its journey from the injection nozzle, its 
position, velocity, temperature and size, in general, will be 
changing according to the interaction between the gas-phase 
solution and the spray properties. In contrast to die Elllerian 
approach for the gas-phase equations, the Lagrangian treat- 
ment for the liquid-phase equations requires interpolation of 
the flow quantities from the Eulerian grid to the particle 
positions and redistribution of the spray source terms from the 
particle positions to the Eulerian grid for the gas-phase. Since 
the time step for die spray equations is usually much smaller 
than that for the flow equations (especially for the present 
MSIP scheme with which, in general, a large CFL number can 
be obtained) and a large number of particle groups are 
required to statistically represent the spray behavior, it is very 
important to keep the numerical efficiency of the present 
implicit scheme from being severely degraded down by the 
coupling (or inter-action) between the gas and liquid-phases. 
A strategy to overcome this problem, at least for steady state 
flows, has been successfully applied to the present study. It 
will be described in this section. 

C.2c(l) Stochastic Process 

As mentioned in the introduction section, there are two types 
of separated flow models for spray computations. One is the 
deterministic separated flow (DSF) model and the other the 
stochastic separated flow (SSF) model. For turbulent flow 
calculations, the DSF model completely neglects the disper- 
sion effect due to the turbulent motion and. therefore, mean 
gas flow quantities are used to evaluate the right-hand-side of 
Eqs. (6) and (7). In the present study, the SSF model is 
applied to account for the dispersion effects on spray charac- 
teristics. In this model, the gas-phase velocity fluctua- 
tions.//, and are generated by randomly sampling a 
Gaussian probability density distribution having a standard 
deviation of (2 k/ 3)* /_ . The instantaneous velocities are then 
used to evaluate Eqs. (6) and (7). A fixed number of samplin- 
gs is conducted for each group of particles. The final spray 
quantities are obtained by averaging the results (source terms, 
trajectories,. . . , etc.) among the total samples. Details of the 
SSF model can be found in Refs. 26 and 33. 
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C.2c(2) Determination of Spray Time Step 

For spray combustion calculations, there are several time 
scales involved in the flow field which can differ by several 
orders of magnitude. The chemical reaction time scale is 
usually very small compared with the rate of evolution of the 
gas flow. Tlus is also true for an evaporating spray. To 
accurately calculate the particle trajectories, size and tem- 
perature. the integration time step has to be small. This is 
especially severe as the droplet becomes smaller and smaller 
toward the end of its lifetime. For an unsteady problem, the 
time step for the entire system will be controlled by the 
smallest time step. For a steady state calculation, however, the 
time step for chemical reaction is usually not a problem in the 
present formulation due to the implicit treatment of the 
chemical source term. Our experience indicates that, for 
single-phase combustion calculations, the same CFL number 
usually can be used for both nonreacting and reacting calcu- 
lations with the present MS1P method. This ensures that the 
convergence properties for reacting flow calculations are not 
degraded using the present numerical algorithm. However, tire 
time step for spray equations still remains small and has to be 
selected (computed) very carefully in order to obtain accurate 
spray results and to maintain stability. The determination of 
the spray time step will be discussed here and how the spray 
interacts widi the gas flow is described in the next section. 

The spray time step at any instant of time along its trajectory 
is determined based on the following time step constraints: (1) 
droplet velocity relaxation time (f r ), (2) droplet life time (f ( ). 
(3) droplet surface temperature constraint time (/,), (4) local 
grid time scale (l g ) and (5) turbulent eddy-droplet interaction 
time (r ( ). The final spray time step (Ar ?pr ) is determined by 
taking the minimum of the above five time steps to ensure the 
accuracy and stability of the spray calculations. A factor 
between 0.1 and 0.5 is further used to multiply the selected 
time step in the current spray calculation. These time steps are 
described as follows. 


In order to ensure the success of the temperature calculation 
using Eq. (13), the exact solution of the infinite conductivity 
model (or uniform temperature model) is used to estimate the 
time step for the present vortex model. The temperature 
equation for the droplet based on this model is 



P/C A 


'( f - rj 




(17) 


Equation (17) has an exact solution (after local linearization) 
in the following form 

= a( 1 - e-*'). 

where A T is the droplet temperature change within the 
integration time step. A, A', and B‘ are defined as 
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Therefore, a time scale can be obtained based on the above 
solution if a desired AT p is specified. This time scale is 
expressed as 



In our spray calculation, AT p * 3 K is specified for those 
particles just leaving the injector, where they experience a 
sudden temperature jump and A T p m 0.5 K for the rest of the 
calculation toward the end of their lifetime. 


Droplet velocity relaxation time (t,.): The local linearized 
droplet equations of motion, Eqs. (6) and (7). have exact 
solutions in terms of the local slip velocity with an expo- 
nential decay form. The time constant for the exact solution 
can be expressed as 


\ 


Pi 

'£ 


V 



Droplet life time (t,): To ensure that the drop size remains 
positive for the practical computational purpose, the droplet 
lifetime at any instant of time is estimated by the following 
equation 


3 ">pPl 

Droplet surface temperature constraint time (t*): When the 
governing equation for the droplet internal temperature 
distribution. Eq. (13). is solved, the temperature solution can 
become completely incorrect due to the use of an inappro- 
priately large time step. This is particularly important for fuel 
with a low boiling temperature (close to room temperature). 


Local grid time scale (lg): A particle can travel across several 
grids and may experience a sudden change of the local gas 
properties if the time step is too large. This not only causes 
inaccuracy in the integration but also increases the difficulties 
of locating the particle positions. Therefore, a time scale. t g . 
is computed to ensure that the particle only moves less than 
one local cell size in one time step. 


Turbulent Eddy-droplet interaction tiuie <tj) : According to 
Shuen et al. (Ref. 42). a particle is assumed to interact with an 
eddy for a time which is the minimum of either the eddy life- 
time or the transit time required for the particle to cross the 
eddy. These times are estimated by assuming that the charac- 
teristic size of an eddy is the dissipation length scale as 


and the eddy lifetime is estimated as 

l = vM) 1/2 - 

The transit time of a particle was found using the linearized 
equation of motion for a particle in a uniform flow 


-T III 


l- ^(tiiT - iT r |) 


(18) 
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where 

8 P l r r 

T = 

3p/dI« ~ « p \ 

and II - 11 ^ . is the relatjye velocity nt the start of the 
interaction. When L e > r | u - n p | , the linearized stopping 

distance of the particle is smaller than the characteristic length 
scale of the eddy and Eq. (18) has no solution. In this case, 
the eddy has captured the particle and the interaction time is 
the eddy lifetime. Therefore 

<i * '<■ L c > T l»* ~ 'V I 
tj = min (r f . r,). if l e < t| IT - n p \ 

C.2c(3) Interaction Between Two Phases 

For the gas-phase equations. Eq. (3). the presence of the spray 
appears in the form of a source term. H ; . This source term 
represents the interchange of the mass, momentum and energy 
between two phases. As the particles are injected into the flow 
domain, their subsequent behavior (positions, velocities, size, 
and temperature) is affected continuously by the neighboring 
gas properties and vice versa. Usually a large number of spray 
particles is desirable to accurately predict the spray behavior. 
However, this requires tremendous computational effort for the 
spray calculation alone. To minimize the computational time 
for the spray, the spray source term. H ( . is not required to be 
updated at every gas-phase iteration (pseudotime time march- 
ing). Usually the spray source term is updated every 10 to 
20 iterations in our spray calculations. When the spray source 
term is updated, each group of particles is integrated either to 
the end of its lifetime or until it leaves the computational 
domain. It should be noted that, for the present steady state 
spray combustion calculation, the spray time step determined 
previously is independent of the pseudotime used for the gas- 
phase equations, which is determined mainly from the obtain- 
able maximum CFL number accord-ing to the local system 
eigenvalues. Therefore, this strategy maintains both the 
efficiency of the flow solver and the accuracy of die spray 
Lagrangian integration. This treatment of the gas-liquid 
interaction is different from those reported by Raju and 
Sirignano (Ref. 17) where time-accurate solutions were their 
primary concern. In the present computation, it is assumed 
that upon impingement with the walls, the droplets evaporate 
completely and assume the local gas flow velocities. The 
interpolation of the gas-phase properties from the Eulerian 
grid to the particle location and the redistribution of the spray 
source term from the particle location to the Eulerian grid is 
applied in a similar way as discussed in Raju and Sirignano 
(Ref. 17). 

CJ NUMERICAL TEST RESULTS 

In this section, results obtained from the ALLSPD algorithm 
with and without spray are presented. These include (1) a 
nonreacting turbulent backward-facing step flow to demon- 
strate the validity of the current turbulence model. (2) a 
single-phase turbulent reverse jet combustion flow to assess 
the present combustion treatment and (3) a spray gas turbine 


combustion flow to qualitatively demonstrate the spray 
calculation for a complex geometry and the interaction 
between the two phases. The spray results shown here 
emphasize the numerical aspect rather than the spray physics. 
The accuracy validation for spray calculations will be consid- 
ered in future calculations. 

C.3a Backward-Facing Step Mow 

The turbulent backward-facing step flow data of Kim et al. 
(Ref. 43) for a two-dimensional channel with an inlet to step 
height ratio of two is selected here to test the validity of the 
k-e turbulence model. A 136 * 100 grid, clustered near the 
step and the top and bottom walls, was used. No chemical 
reactions were included in the calculation. 

The particle traces of the flow is shown in Fig. 8. The 
experimental reattachment length given by Kim et al. is 7.1 H, 
where H is the step height. The predicted value in our calcula- 
tion is about 6.1 H. which represents a 14 percent under- 
prediction. Figure 9 shows the mean velocity profiles at 
various axial locations. The agreement is quite good in all 
locations except near the reattachment point. The profiles for 
the turbulent kinetic energy and turbulent shear stress are 
given in Fig. 10. Good agreement is observed in these com- 
parisons. Wall pressure coefficients along the step-side and 
opposite walls are shown in Fig. 1 1 . The underprediction of 
reattachment length is also evident from this figure. The pre- 
dicted pressure recovery downstream of the reattachment point 
is in excellent agreement with measured pressure levels. The 
convergence history for this calculation is illustrated in 
Fig. 12. The convergence property for the calculation is 
satisfactory. 

C.3b Reverse Jet Combustion How 

The flow configuration is a 51 mm l.D. (inner diam) times 
457 mm cylindrical chamber containing a reverse jet flame 
holder which issues from a 1.32 mm l.D. (6.35 mm O.D.) 
tube. The jet is coincident with the chamber axis and located 
80 mm upstream from the chamber exit. Both the main and jet 
flows are stoichiometrically premixed propane and air at a 
temperature of 300 K. with a mean velocity of 7.5 m/s for the 
main stream and 135 m/s for the jet. A complete description 
of the flow system is available in McDannel et al. (Ref. 44). 
A 1 17*61 grid (half domain) clustered near the jet tube was 
used. Five species (CjHg, 0 2 , N 2 , C0 2 . and H 2 0) were 
considered in this calculation and the single-step global 
reaction chemistry model reported in Westbrook and Dryer 
(Ref. 45) was used for combustion. Figures 13. 14. and 15 
show the particle traces, velocity vectors (colored by tem- 
perature) and temperature contours of the reacting flow, 
respectively. The experimentally measured temperature con- 
tours <directl> taken from McDannel et al. (Ref. 44)) are also 
presented in Fig. 15 for comparison. The flow is clearly seen 
to consist of two distinct regions - the recirculation zone and 
the wake. The incoming flow is ignited by the hot combustion 
gas in the recirculation zone and further combustion takes 
place in tine wake. For the test conditions considered here the 
reverse jet serves as a very effective flame holding device. 
The results in Fig. 15 indicate that the predicted temperatures 
are higher than the measured values. This is mainly attributed 
to the over-simplified chemistry model used in the present 
calculation. As reported by McDannel et al. (Ref. 44). there 
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was a significant amount of carbon monoxide (CO) observed 
in the combustion products, which lowered the combustion 
temperature. Also not included in the calculation is the radia- 
tion heat loss which will lower the flow temperature further. 
Ina addition, the thermocouple's used in the measurements 
were not corrected for radiation loss, indicating that the actual 
gas temperatures would be somewhat higher then measured. 

C.3c Gas Turbine Spray Combustion Flow 

A simplified model of General Electric's EEE (Energy 
Efficiency Engine) combustor is considered. Figure 16 shows 
the clustered 8 1 *65 grid for this calculation and the general 
engine grid layout for the entire engine (Ref. 46). It is an 
annular combustor with a dual can combustor dome. The 
cooling air through the internal walls (two combustor domes) 
was omitted for simplicity. Also, the swirling effect was not 
included in the calculation. The present calculation does not 
completely simulate the typical gas turbine combustion 
characteristics, in which swirling and cooling are the two 
important ingredients. We would like to emphasize, in the 
present study, that the focus is to demonstrate the effective- 
ness of the interaction between the two phases based on the 
present spray solution procedure. Further studies with this 
algorithm will be conducted both on the spray accuracy 
evaluation and on the detailed swirling and cooling computa- 
tion for gas turbine combustor configurations. 

For the present spray combustion calculation, the spray 
injectors were located close to the inlet of the combustor 
domes. One hundred spray groups with 10 random samples for 
each groups were used, which can be seen in Fig. 17. The 
case studied here has a flow Reynolds number. Re * 1 .05* 10*. 
where the Reynolds number is based on the inlet maximum 
velocity and combustor inlet height. The inlet air temperature 
is 900 K and the pressure is I atm. The liquid n-pentane fuel 
was used and five species ( CjH (? 0 2 , N 2 , C0 2 , and H 2 0 ) 
were considered in this calculation. Again, the single-step 
global reaction chemistry model reported in Westbrook and 
Dryer (Ref. 45) was used for combustion. The fuel/air ratio is 

0.02 (total fuel injected/total incoming air at the inlet, includ- 
ing bypass air). The injection velocity of the liquid fuel at the 
exit of the fuel nozzle is 20 m/s and die temperature is 290 K. 
The liquid fuel was assumed to be fully atomized with the 
initial diameters ranging from 20 to 100 pm. The liquid fuel 
was injected into the gas flow after the gas flow had been 
iterated to reach a nearly steady state solution. Upon the 
injection of die fuel, the cool fuel was suddenly exposed to a 
hot environment and the interaction between the two phases 
took place in terms of die interchange of mass, momentum 
and energy. An ignition source was placed downstream of die 
injector to ignite die "buming" of the fuel-air mixture. The 
ignitor was turned on right after the initiation of the spray and 
was turned off when the temperature in any of the ignition 
computational cells reached 1600 K. The spray source term. 
H,. in Eq. (3) was updated (computed) every 20 gas-phase 
iterations. The converged results are presented in the following 
figures. 

Figure 10 shows the particle trajectories. The dots in die 
figure are not scaled to indicate the relative size of the actual 
liquid particles, although smaller particles evaporate much 
faster than bigger ones. The computed velocity vectors colored 


widi temperanire are presented in Fig. 18. where a converged 
solution of the cold flow (nonspray/noncomhustion) with the 
same flow condition is also included for comparison. The 
combustion zones are confined near the inner walls. This is 
mainly due to the exclusion of the swirlers at the inlet of the 
combustor domes in the present calculation. Without swirl, the 
degree of the fuel-air mixing is relatively poor. As may be 
noted, there is no flame holding device for the present 
calculation. The recirculation zones near the inner walls are 
the only devices to hold the flame. Figure 19 shows the Sauter 
mean diameter. D J2 . along die averaged trajectory for both the 
lower and upper domes. The increase of the D 32 right after the 
injection indicates the rapid evaporation of the smaller 
particles. After reaching a peak, both mean diameters gradu- 
ally decrease as more and more particles evaporate. To show 
the convergence properties for both nonreacting and spray 
combustion cases, the convergence histories for both cases are 
shown in Fig. 20. As can be seen in this figure, the spray was 
initiated after the gas-phase solution had been iterated 1000 
times. A sudden disturbance from the spray injection causes 
the L2 norm residual to jump to a level higher than die initial 
gas-phase residual. As more and more spray particles evapo- 
rate and undergo combustion, the interaction between two 
phases can be clearly identified in this convergence pattern. 
Since in the present calculation, the spray source term was 
updated every 20 gas-phase iterations, a small residual spike 
along the convergence history can be seen very clearly. These 
small spikes persist toward the end of the present computation. 

D. CONCLUDING REMARKS 

Much work is required to develop a CFD code to the level 
necessary for a designer to be able to use it with confidence. 
We at NASA are proceeding with a commitment to do this. A 
three-dimensional version of our ALLSPD will be developed 
shortly and extensive validation against experimental data will 
be performed. At the same time we are exploring the use of 
advanced modeling to better represent the flow physics and 
chemistry of turbulent combustion. For example, modeling of 
dense sprays and sprays at super-critical conditions will be 
investigated. Modeling of the turbulence and chemistry inter- 
action will be explored, using techniques such as PDF 
methods. The use of reduced hydrocarbon chemical kinetics 
models which represent the actual combustion processes are 
also necessary. And finally, an improved model to represent 
potentially high levels of radiation heat transfer is also needed. 

The state of art for numerics also continues to advance at a 
rapid pace. Our code will be modified in the future as oppor- 
tunities arise in areas such as unstructured grids and massively 
parallel computing. This ongoing effort to produce a modem 
chemical reacting flow CFD code holds promise to provide a 
powerful design tool for the industry's use in the analyses of 
the next generation of gas turbine engine combustor concepts. 
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Figure 2. — Optical configuration of the phase/Doppler 
particle analyzer. 
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Figure 3. — Gas phase mean velocity at 5 mm down- 
stream. 
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Figure 4. — Gas phase fluctuating velocity at 5 mm down- 
stream. 
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Figure 6. — Drop fluctuating velocity at 5 mm downstream. 


Figure 5. — Drop mean velocity at 5 mm downstream. 
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Figure 7. — Drop number flux at 5 mm downstream. 



Figure 8. — Particle traces for turbulent backward-facing step flow. 
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Figure 9. — Numerical predictions of the mean velocity profiles at selected axial locations for the turbulent 
backward-facing step flow. 
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o Experiment of Kim et al., 1 980 
Computations 
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(a) Turbulent kinetic energy. 
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(b) Turbulent shear stress. 


Figure 10. — Comparisons of (a) turbulent kinetic energy and (b) turbulent shear stress at some selected axial 
locations for the turbulent backward-facing step flow. 
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o Experiment of Kim et al., 1 980 
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(b) Opposite wall. 


Figure 11 . — Comparison of the wall pressure coefficients 
for the turbulent backward-facing step flow. 



Figure 1 3. — Particle traces for the reverse jet combustion 
flow. 



Figure 1 4. — Velocity vectors (colored by temperature) for 
the reverse jet combustion flow. 



Figure 12. — Convergence history for the turbulent backward- 
facing step flow. 



(a) Experimental data (from McDannel et al.). 



Figure 1 5. — Comparison of the temperature predictions 
with the experimental data for the reverse jet combustion 
flow. 
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Figure 16. — Center-plane cut for GE's EEE gas turbine 
engine and the grid for the 2-D combustor. 


Figure 1 7. — Liquid fuel particle trajectories for the 2-D 
gas turbine spray combustion flow. 
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(a) Non-combustion flow (cold flow, no spray). 


(b) Spray combustion flow. 




Figure 1 8. — Velocity vectors (colored by temperature) for (a) non-combustion and (b) spray combustion cases for the 2-D 
gas turbine spray combustion flow. 
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Figure 1 9. — Sauter mean diameter distribution for upper 
and lower combustor domes along the averaged tra- 
jectory for the 2-D gas turbine spray combustion flow. 



Figure 20. — Convergence history for both non-combus- 
tion and spray combustion cases for the 2-D gas tur- 
bine spray combustion flow. 
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